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We perform a variational Gutzwiller calculation to study the ground state of the repulsive SU(3) 
Hubbard model on the Bethe lattice with infinite coordination number. We construct a ground-state 
phase diagram focusing on phases with a two-sublattice structure and find five relevant phases: (1) 
a paramagnet, (2) a completely polarized ferromagnet, (3) a two-component antiferromagnet where 
the third component is depleted, (4) a two-component antiferromagnet with a metallic third com- 
ponent (an "orbital selective" Mott insulator), and (5) a density- wave state where two components 
occupy dominantly one sublattice and the last component the other one. First-order transitions 
between these phases lead to phase separation. A comparison of the SU(3) Hubbard model to the 
better-known SU(2) model shows that the effects of doping are completely different in the two cases. 

PACS numbers: 67.85.-d,71.10.Fd,64.75.Gh 
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I. INTRODUCTION 

The recent success in loading an ultracold fcrmionic 
cloud in an optical lattice with a total of eight 
fermion components and SU(2) x SU(6) global sym- 
metry using alkaline-earth- metal atoms [l| demonstrates 
that fermionic SU(-/V) Hamiltonians can be realized ex- 
perimentally. A simple way to create SU(iV) symmetric 
models is to use the fact that neither the kinetic nor the 
interaction energy of atoms with closed electronic shells 
depend on the orientation of the nuclear spin. For a given 
nuclear spin / one can therefore realize SU(A'^)-symmetric 
models with < 2/ -1-1 by selectively occupying A^ of the 
2/ -1-1 spin states [2|. For deep enough optical lattices the 
system is naturally described by a Hubbard model with 
a good accuracy Q 13] • Generalizations of the Heisenberg 
and Kondo Hamiltonians have also been proposed Q to 
describe the system in other parameter regimes. While 
the N = 2 Hubbard model has been studied extensively 
(see and references therein) and many results were 
obtained in large- A^ expansions @ (mainly focusing on 
even values of A^ and the vicinity of half filling), much 
less is known for the SU(A^) case (far) away from half 
filling and if A^ is odd. Experimentally, the first indica- 
tions of the Mott transition in the SU(2) Hubbard model 
have recently been observed using fcrmionic atoms in op- 
tical lattices 0, Q but lower temperatures are needed to 
stabilize, for example, an antiferromagnetic phase. 

In this work we shall focus on the ground state of 
A^ = 3 fermionic components (colors) in a lattice with 
isotropic, local interactions. The case of attractive inter- 
actions have already drawn considerable attention [ol-fl^ 
due to similarities to quarks [l3| . however, many-body 
losses (e.g., [lj| for ^Li) and high temperatures make 
experimental observation of the corresponding color su- 
perfluid phase and the liquid of the three-body bound 
states ("baryonic" phase) difficult. In the case of the 
repulsive interactions, Szirmai, Legeza, and Solyom [l^ 
used DMRG and bosonization to study the phase dia- 
gram in one dimension, D = 1. Already some time ago 



Honerkamp and Hofstctter [l6| used functional renormal- 
ization group (fRG) in the D = 2 square lattice close 
to half filling to investigate the interplay of spin- and 
charge-density waves and the possibility of staggered flux 
phases for larger A^. More recently, Gorelik and Bliimer 
|17| performed paramagnetic dynamical mean field the- 
ory (DMFT) calculations to discuss the Mott phases. 
Finally, Miyatake, Inaba, and Suga [l^ studied two- 
sublattice ordering and Mott transitions at half filling as 
a function of the anisotropy in the interactions between 
the components using DMFT. Nevertheless, a compre- 
hensive study of the phase diagram for arbitrary fillings 
is still missing. We apply a non-perturbative method, 
Gutzwiller ansatz with Gutzwiller approximation, to in- 
vestigate the T = phase diagram which could serve as a 
starting point for future study. The Gutzwiller approach 
is a method which captures the physics both at weak 
and strong coupling and has strongly influenced our un- 



22 . We concentrate 



derstanding of correlated matter \li 
on phases with two-sublattice ordering. 

We will approximate the ground state of the Hamilto- 
nian 



(ij)\a 



(1) 



where a — \, 2, 3 denotes the three fermionic compo- 
nents, i denotes sites in a bipartite lattice with z nearest 
neighbors, t* is a normalized nearest-neighbor hopping 
amplitude, and U > is a local repulsive interaction. 
We shall calculate ground state properties as a function 
of U and the filling p, where p = 1 corresponds to total 
filling, that is, three fermions per site. 

The Hamiltonian in Eq. ([T]) has a global SU(3) sym- 
metry due to the isotropic interaction strengths between 
the components. Our main goal is to study how this sym- 
metry breaks in the ground state due to the competition 
between kinetic and interaction energy. The scaling of 
the hopping ensures that the kinetic energy of Eq. ([T]) 
per lattice site becomes independent of z in the limit 



2 



2 — >■ (X). As we will see, the Gutzwiller approximation, 
used below, is exact in this limit [l^ . 

The Hamiltonian ([T]) can easily be generalized to 
SU(-/V) global symmetry using N components instead of 
three. We will mainly study the TV = 3 case, but consider 
arbitrary N in some cases. 

It is easy to see that the model has a particle-hole sym- 
metry as the canonical transformation Cia — >■ (~1)*c|q, 
leaves the Hamiltonian invariant up to a constant term. 
Thus the phase diagrams for p < 1/2 and p > 1/2 can 
be mapped onto each other. 



II. GUTZWILLER ANSATZ 

We approximate the ground state of the Hamiltonian 
in Eq. ([1]) by a Gutzwiller wavefunction. 



(2) 



The unprojected wavefunction |^'o) is a normalized 
Slater determinant. For simplicity, we assume that 



(3) 



that is, we consider only collincar order where the mag- 
netization matrices {'i'o\cf^Cip\'i'o) on different sites com- 
mute. 

We shall use the notation 



>o|^i^a|^'o) 



(4) 



The operator pj projects on the local configuration 
/ e {0, 1, 2, 3, 3 = 12, 2 = 13, T EEE 23, t = 123}, satisfying 
PiPi' = Sii'Pi- The Gutzwiller parameter therefore 
changes the relative amplitude of the configuration / at 
site i. It is useful to define occupation operators n/, with 
■hinp = niup. The two operator sets are related by 

^p,, ;j5,= ^(-l)l^'H^ln,,. (5) 

I'DI I'DI 

The Gutzwiller expectation value of an operator O is 
defined as 



o = 



{G\d\G) 



(6) 



(G|G) ■ 

We shall calculate the variational energy, defined by 

= (7) 
and minimize it with respect to |^'o) and A/(i). 

III. GUTZWILLER APPROXIMATION 

Unfortunately, the exact analytic evaluation of the 
Gutzwiller expectation values in Eq. ([7]) is in general not 



possible. Originally, Gutzwiller [20[ proposed a mean- 
field-like approximation to calculate the variational en- 
ergy. Later it was shown in the N = 2 case that in 
the limit of infinite coordination numbers, z — >■ oo, the 
Gutzwiller approximation is exact [l^, . Here we shall 
outline an approach to calculate the variational energy 
in Gutzwiller approximation in the A'^ = 3 case based 
on functional integrals. Technical details can be found 
in the Appendix. An alternative method to derive the 
variational energy within the Gutzwiller approximation 
has been described for multiband models by Biinemann, 
Gebhard, and Weber in Ref. [13 . 

The main idea is to express all Gutzwiller expectation 
values as expectation values in a static auxiliary field 
theory, as in Ref. Q . The simplest way to construct this 
field theory is to start from the norm of the Gutzwiller 
wave function. We observe that after normal ordering 
and then using Wick's theorem we can always write 



(G|G> 



n 



, /'C/ 



,(8) 



where [Go]ia-ji3 = (^o|c]^CiQ,|^'o) is the bare equal-time 
propagator of the unprojected wave function, = Cia 

is a time-independent Grassmann field, and the occu- 



pations, fiij 



n 



are expressed in terms of 



the Grassmann variables Cia- Using the Grassman alge- 
bra, we can exponentiate the terms in the product, and 
rewrite Eq. ([5]) as a "partition function" , 

^aux[*^*] = 'fH-Go')^ + Y.ui{i)nu, (9) 

i-j 

where ui{i) are functions of the Xj{i). It can be seen 
that for any operator, (G|0|G) can be expressed as some 
functional integral with the same static action S'aux- As a 
consequence, Gutzwiller expectation values become cer- 
tain expectation values in an interacting field theory. We 
can calculate all such expectation values in the limit of 
z — >■ oo using the cavity method [1^ described in the 
Appendix. Here we use that to obtain the Gutzwiller ap- 
proximation for the static field theory S'aux, one follows 
the same steps which have to be taken to derive dynam- 
ical mean field theory (DMFT) [25| in the large z limit. 
As there is no time dependence, the relevant Grassmann 
integrations become finite dimensional and can be cal- 
culated analytically (see the Appendix). Note however, 
that the coupling parameters of Saux remain variational 
parameters and one still has to minimize the resulting 
energy. 

In the limit of z — > oo. the proper self-energy matrix 



Gn 



-^'^''1')^^ corresponding to the action S'a 



becomes site diagonal 



as the arguments 
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of Ref. [l9[ can be generalized to include higher order 
local interaction terms. It is important to note that in 
contrast to DMFT, the self energy is static, that is, it 
does not depend on frequency. Therefore all self-energy 
effects can be absorbed into Go. Thus, by restricting the 
variational parameters Xg}{i) and Xa(i) such that 



(G|G> = 1 and 



(^'olnial^o) 



(V«), 



(10) 
(11) 



one can set the self energy to zero (Sic = 0). 

Note that in Eq. ([TT|) we simply require that the Gutz- 
willer projection has to leave the densities invariant. As 
a consequence, the parameters Xa{i), a — 1, 2, 3, are 
simply replaced by the order parameters of |5'o) as vari- 
ational parameters. Furthermore, it is also possible to 
replace the two- and three-body Gutzwiller parameters, 
Xa{i) and At(z), by the physical occupation numbers, 
dia = n-ia and ti = na as variational parameters (see 
the Appendix for details). 

After some algebra, the variational energy can be cast 
to the form 



E-i, 



H.c. 



qta Qja (fo I cLcja I *o) 

+ f/^rf.a, (12) 



unprojected wave function |^'o). As in Ref. 12g, one can 
perform the variation with respect to (^ol analytically. 
This leads to an effective mean-field Schrodinger equation 



i^ol*o) = i?sp|*0> 



(16) 

(ij);Q! ia 

with I^I'o) being the ground state of Hq. Therefore the 
unprojected wave function |^o) is a Slater determinant 
as required by our approach. For a specific set of {jt-^q,} 
(these also define in our case the order parameters in the 
unprojected state), one could in principle construct the 
unprojected wave function. However, we do not need 
1 4*0) explicitly, only the corresponding ground-state ex- 
pectation values. For this reason, we shall perform cal- 
culations on the Bethe lattice with infinite coordination 
number, where the (ground state) Green's function cor- 
responding to Hq can be calculated analytically. Another 
advantage of doing so is that on the homogeneous Bethe 
lattice, the density of states (DOS) is semieUiptic. This, 
on one hand, is a good approximation to the density 
of states of a d = 3 simple cubic lattice with nearest- 
neighbor hopping, while on the other hand, it allows us 
to compare our results quantitatively to DMFT studies, 
which used the same bare DOS (see Refs. [l3l and [l8|). 



where the renormalization factors are given by (shown 
here for a = 1, other components are related by permu- 
tations) 



v/<(l-0 



(13) 



Here, = {G\pii\G) / {G\G) are functions of the physical 
occupancies nu [see Eq. ([5])]. Note that for any physical 
solution, has to be real, and thus there is a number 
of constraints on the variational space corresponding to 



P^I > 0. 



(14) 



Equation (|T2|) is consistent with the general form of 
the variational energy in Gutzwiller approximation ob- 
tained by using other techniques (see Refs. [13, H^). The 
physical interpretation of the renormalization factors qia 
is that they describe how the noninteracting band gets 
renormalized by the interactions. 

Finally, for a fixed set of the bare occupations {n^^}, 
the variational problem can be written in the form 



min [S„ -£;,p((*o|*o) - 1) 



(15) 



Here Esp and Xia are Lagrange multipliers and is a 
function of j^'o) and the double and triple occupancies, 
dia , ti- We observe that this expression is quadratic in the 



TWO-SUBLATTICE ORDER ON THE 
BETHE LATTICE 



IV. 



The variational problem in Eq. (jisp can be applied to 
describe collinear magnetic structures with an arbitrary 
large magnetic unit cell. However, the increasing num- 
bers of the variational parameters and, especially, the 
constraints [Eq. (jl4| . makes evaluation and minimiza- 
tion more and more difficult as the size of this unit cell 
grows. Thus we shall concentrate on structures with two- 
sublattice order. The relevance of other phases is dis- 
cussed in the concluding section. 

The on-site components of the ground-state Green's 
function on the Bethe lattice with two inequivalent sub- 
lattice in the limit z — > cx3 are given by [27| 



Gaii,i;uj) 



(17) 



where W = 2i*, lia = (w + Xia)/qfa- In the case = 
and qia ~ 1 one simply recovers the semielliptic density 
of states from the jump of the Green's function at the 
branch cut, 



(18) 



It is useful to measure energies from the "chemical poten- 
tials" fia = —{Xia + Xi+ia)/'2, and introduce the energy 
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difference between A and B sites, ha = |AiQ, — Ai+ia|/2, 
and the renormalized bandwidth, Wa = WqiaQi+ia- In 
terms of these parameters, the two subbands which arc 
given by the branch cuts of the Green's function in 
Eq. p7|) he in the intervals fia 

and Ha + ha < LU < fla ^ 



Vhl + Wl < U < 



Using the Green's function in Eq. p7|) it is possible 
to express the local densities as (^'ol^ml^'o) = Pa + 
(— l)*mQ,/2. We can get relatively simple expressions for 
the homogeneous and staggered parts, namely 



Pa 



de Do{e) 



I I epa L |eFap 1 . £Fa .,„^ 

— \ I — — — arcsm (19) 



and 



— "m "-i+la 



de Doie) 



2Aa 



-w 



TT W 



1 + ^ {Fiea,ka) - E{ea, ka)) ,(20) 



M/2 



where 



arccos(| — 



and 



^\)^ka ^ 1/^1 + /^l/W\ 

the functions F and E are incomplete elliptic inte- 
grals of the first and second kind (defined according 
to Abramowitz and Stegun [1^), respectively, while the 
roles of the parameters tpa and are discussed below. 

The variational energy for two-sublattice long-range 
order can finally be written as 



Ev^^i QAaQBaKa + U 



dAa 



(21) 



where the mean-field kinetic energy can also be expressed 
in a compact form. 



de Doie) 



-w 

2W 

2W 
"3^ 



-(1 



A2 



A2 

E{ea,ka) 



. / 1 ^Fa \^Fa\ / -^pa , ~a /r,r,\ 



Equations (|19p . ((20)) . and ((22|) generalize the expressions 
corresponding to the two-component antiferromagnetic 
ansatz in Ref. [l^ to the three-component case on the 
Bcthe lattice. 

The variational energy for a given p and U and a 
given set of variational parameters, pa (with a constraint 
J2aPa — 3p), Aq,, dia and ti (i = A,B), can be calcu- 
lated as follows: First, one has to solve Eq. ([T^ for each 



a = 1, 2, 3 for epa- Then calculate ttIq and Ka from Aq, 
using Eqs. (PH]) and (f^ . respectively. Last, the occupa- 
tions 77°^ = Pa + (— 1)*777q/2, dia, and determine the 
renormalization factors qta- In principle, one could use 
the nia instead of the A^ as variational parameters, but 
then one has to invert Eq. ([20]) . Finally, the energy has 
to be minimized taking also the constraints into account. 
Since the minimization of Eq. (j2ip with respect to the 
variational parameters cannot be performed analytically, 
in the next section we discuss the results of the numerical 
optimization. 



V. RESULTS 

To explore the phase diagram in the general case with 
two-sublattice symmetry, we added a quadratic penalty 
function ~ (^^ -^^ ^ Q{—Pii)Pii)'^ to the variational en- 
ergy and minimized it using different stochastic methods 
which gave consistent results. We found five phases which 
we discuss in the subsections below. The sketches of the 
occupations of the corresponding phases on the two 
sublattices are shown in Fig. [T] Note that the minima 
are often found at the boundary of the variational space 
defined by the constraints. 
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FIG. 1: (color online) Examples of the sublattice occupations 
of the phases discussed in the text. The different colors (blue, 
green, red) represent the different components (a = 1, 2, 3, 
from left to right), the heights are proportional to the occu- 
pations. 



To check whether the correct global (rather than lo- 
cal) minima have been found, we also performed a brute- 
force random Monte Carlo search (MC) in a restricted 
variational space where we simply threw away configura- 
tions which violated the constraints. We assumed that 
Pi = P2 = P + mo/6, p3 = p - mo/3, 7711 = -777,2 = mq, 

TO3 = 0, dAl = dB2 = di,dA2 = dsi = d2,dA3 = dB3 = 

c?3 and tA ^ tB ^ where A and B refer to the two 
sublattices. Despite the significant reduction of the vari- 
ational space, four of the five phases can be described 
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using this parametrization. In Fig. [2] a typical result of 
such a search is shown. 

For a precise calculation of the location of the phase 
transitions we performed in a third step variational cal- 
culations for each phase separately using the appropriate 
subset of variational parameters. 

In the following we will describe each of the five rele- 
vant phases which arc also sketched in Fig. [TJ 
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FIG. 2: (color online) Raw output of parameters mo — pi + 
p2 — 2p3 and mq = mi = — m2 after a MC search for the 
minimum for different values of p for U /W = 5. For each 
value of p we used 2^" random points in variational space. 
We see three regions with distinct behavior, corresponding to 
the PM, AF2, and AFMM phases. 



A. Paramagnetic state (PM) 

At low enough values of p or [/, no ordering is expected 
and the ground state has to be a (correlated) paramag- 
net, defined by toq = TOq, = 0, dia = d, U = t. Such a 
paramagnetic Gutzwiller ansatz can also be extended to 
the SU(A^) case with N components, giving in Gutzwiller 
approximation 



E:, 



PM 



(23) 



Here Ne is the total kinetic energy of the non-interacting 
Fermi sea with filling p. In a simple, low density approx- 
imation (by neglecting high-order occupancies, n\j\ = 
for |/| > 2), we find a relatively simple expression for the 
renormalization factor 



Q 



PM 



{N -l)d 



^l'Np+ (^^^ d+{N-l)Vd . 



(24) 



Note that this expression, which neglects triple and 
higher occupations, is exact for N = 2. It also gives 



a good upper bound for the ground-state energy for low 
densities, p < 2/iV, becoming even better as U increases. 
The advantage of using this approximation is that we can 
obtain analytic expressions for general N. 

In particular, we see that at the commensurate fill- 
ing p = 1 /N the variational energy becomes a quadratic 
function of d, ~ d{d—do) with d > 0. The sign change 
of do ^ Ubr ~ U signals the Mott transition (found by 
Brinkmann and Rice [2l| for N — 2). The Mott phase 
appears for 



[/ > C/br = 2 




TV 



2(iV- 1) 



Ne. 



(25) 



For N — 2, this reduces to the well-known result [2l|, |29|, 
while for N = 3 it gives UsRiN = 3) w 3.975VF, which 
matches (up to a relative error of lO"'') the numerically 
exact value obtained from Eq. (|23|) when one restores the 
triple occupancy i as a variational parameter, using 



iPM/ 



(iV = 3) = V(l -Sp + Sd-t)ip-2d + 1) 



+2^/{p-2d + t){d-t) 



+Vid^t /b(i-p)]- 



(26) 



Particle-hole symmetry implies that there is another 
Mott phase at p = 1 - 1/3 = 2/3. 

The results of the numerical minimization of Eq. (|23p 
with Eq. (j26p with respect to d and t are displayed in 
Fig. [3] which shows the filling as a function of the chemi- 
cal potential. The Mott insulator found for p = 1/3 (and 
for p = 2/3) is characterized by a vanishing compressibil- 
ity, dp/dp = 0. It is instructive to compare these results 
to the DMFT results of Gorelik and Bliimer [13 who 
found a Mott insulator for U > 2.75W. The quantitative 
discrepancy to our T ~ result can, however, be traced 
back to the finite temperature, T/W = 1/40, used in 
the DMFT calculation. Using the Kotliar - Ruckenstein 
slave-boson mean field theory [s^l, which is a natural 
generalization of the Gutzwiller approximation to finite 
temperatures, for the N = S paramagnetic case we ob- 
tain (not shown) at T/W = 1/40 a critical U « 2.83H^, 
close to the DMFT value. Despite this good agreement, 
one should, however, keep in mind that DMFT and the 
Gutzwiller calculation provide very different scenarios of 
how precisely the Mott transition occurs [25j . In reality, 
however, the T — Mott transition is masked by various 
ordered phases discussed below. 



B. Completely polarized ferromagnet (FMl) 

For the two-component Hubbard model, it has been 
proven by Nagaoka in Ref. [3l[ that for a class of lat- 
tices the exact ground state for U /W — > oo and a half- 
filled system with a single hole is a completely polar- 
ized ferromagnet. The physical argument is that all spin 
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FIG. 3: (color online) The filling p as a function of the chem- 
ical potential /i = dE^^ /3dp in the paramagnetic Gutzwiller 
calculation, for U/W = 0, 1, 2, 3, 4, 5. We see clear Mott 
plateaus at p = 1/3 and p = 2/3 for U/W > 3.975. Inset: the 
double occupancy d as a function of the filling p. 
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FIG. 4: (color online) Stability regions of homogeneous phases 
with two-sublattice symmetry based on the Gutzwiller calcu- 
lation without taking into account that phase separation can 
occur (see Fig. [5] for comparison). The inset shows the tiny 
region to the left of the dashed line where one obtains a par- 
tially polarized ferromagnet with a polarization close to (but 
not exactly) 100%. All solid lines are first-order transitions. 



configurations are exactly degenerate in the exactly half- 
filled system for [/ = 00 and the kinetic energy is fully 
quenched. Therefore one has to ask the question which 
spin configuration allows for a maximal gain of kinetic 
energy for an additional particle or hole. For N ~ 2 and 
p = 1/2 the maximal kinetic energy — is gained only 
for a ferromagnetic arrangement of spins. For = 3 and 
hole doping of a p = 1/3 state, that is, for the removal 
of a particle, the same physical argument applies but it 
is not valid when an extra particle is added instead. For 
example, a 'red' particle added to a 'blue-green' antiferro- 
magnetic configuration also gains — in kinetic energy. 



In this case, higher order terms of order W"^ /U which 
favor antiferromagnetism will suppress ferromagnetism. 

Accordingly, we find within our Gutzwiller calcula- 
tion that for very large U a fully polarized ferromag- 
netic state exists for fillings 1/4 < p < 1/3 (see Fig. |4]). 
Only in a small window of parameters, for example, 
0.25 < p < 0.255 for U/W = 00, (see inset of Fig. S]) 
we find a partially polarized ferromagnetic state with a 
very large value of the polarization which jumps to zero 
at the first-order transition to the paramagnetic state. 
Note that this result is sensitive to details of the density 
of state. 

The discussion can be generalized to arbitrary N as- 
suming a first-order transition from the paramagnet to 
the fully polarized ferromagnet. The energy of the ferro- 
magnetic state, FMl, is simply given by 



E 



FMl 



(p) = Jde Doie)e, 

-w 



(27) 



which has to be compared to the paramagnetic energy 
discussed above. For U ^ 00, when there are only singly 
occupied or empty sites, the energy of the paramagnetic 
state is obtained by setting d = t = in Eqns. and 



£;™(p) 



1 - Np 



N de Doie)e. 



(28) 



-w 



Both E™'^ and E^^ vanish for p=l/N when the para- 
magnetic solution describes a Mott insulator and the fer- 
romagnetic one a band insulator. Furthermore, for any 
symmetric density of states, the two energies also coin- 
cide for p = 1/(A^ + !)• In this case the filling of the 
ferromagnetic band is N/{N + 1) = 1 — 1/{N + 1) and 
thus the integrals in Eqs. ([27)1 and (pS]) . where the band 
filling is 1/(A^ + 1), coincide. Also the renormalization 
factor in Eq. ((28|) is 1/A^ and cancels with the factor 
N arising from the N bands. Therefore one finds that 
within the Gutzwiller approximation of the SU(A^) Hub- 
bard model the fully polarized ferromagnetic state has a 
lower energy than the paramagnet for 



N 



1 1 
< p < — 



(29) 



in the limit of ?7 — > 00. 

To check for the possibility of a partially polarized 
state, we also calculate for J7 — >■ 00 the energy differ- 
ence, -Bflip, when a single particle has changed its color 
compared to the fully polarized state. 



Emp = -(1 



„FM 



)W- 



1-p 



FM 



de Do(e)e, 



-w 



(30) 
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with p™ = deDo{e) = N p. Here the first two 

terms describe the change of kinetic energy of the flipped 
particle and the third term reflects that the kinetic energy 
of all the other particles gets renormalized by the color 
flip. For e^j^/W > 0.4317, Efnp is positive, implying that 
the completely polarized state is stable. For = 3 and 
p < 0.255, however, spin flips are energetically favored, 
as discussed above. For > 4 on the other hand, the 
completely polarized state remains stable in the whole 
interval described by Eq. (|29)) iov U/W oo. 

C. Two-component color antiferromagnet metal 
and insulator (AF2) 

For p < 1/3 filling and moderately strong U, the lowest 
energy state according to the numerical minimization can 
be characterized by mi = — m2 7^ 0, 7713 = 0, ma = 'Sp, 
^ 0, implying pa = and di ~ d2 ~ t — (see 
Fig. 2]). This describes (see Fig. [T]) an antiferromagnetic 
metal where only two of the three components occur. 
With these parameters, formulas are analogous to the 
SU(2) antiferromagnetic state discussed in Refs. [l9l.[2^. 
with a different density of states, however. At p = 1/3 
one obtains a two-component antiferromagnetic insula- 
tor, which has a lower energy than the three-component 
PM at p = 1/3 if t/ > [/ncoI ~ 2.107 VK. Below we will 
discuss that the two-component antiferromagnetic metal 
is never realized as a lower energy is obtained by phase 
separation into the antiferromagnetic Mott insulator and 
either a paramagnetic state or ferromagnetic state with 
lower density. 



D. Three component color antiferromagnetic metal 
(AFMM) 

For 1/3 < p < 1/2 and for sufficiently strong U, the 
dominant state has parameters mi = — m2 7^ 0, ~ 
0, mo = 3 — 6/9. Thus the first two components have 
exactly commensurate filling, pi = p2 = 1/2, forming an 
antiferromagnetic insulator, while the third component 
remains metallic. This phase could thus be called an 
"orbital selective" Mott insulator. At p = 1/2, DMFT 
results are also available for this phase from the work 
of Miyatakc, Inaba, and Suga For C/ = 2.51^, p = 

1/2, the energy and the order parameter mi of the two 
methods agree within a few percent. The AFMM state 
was also found by Honerkamp and Hofstetter [l^ close 
to p = 1/2 on a square lattice using fRG. 

Both AF2 and AFMM phases can be visualized as be- 
ing simultaneously ferromagnetic and antiferromagnetic. 
Here "ferromagnetic" order means that two components 
are equally populated (pi = p2 in our parametrization) 
while the third one is different, breaking SU(3) symmetry. 
The remaining SU(2) symmetry within the first two com- 
ponents is broken by their antiferromagnetic ordering. 
Note, however, that this type of SU(3) symmetry break- 



ing should not be confused with ferrimagnetism, that is, 
simultaneous ferromagnetic and antiferromagnetic order 
in the SU(2) case. For SU(2) symmetry, the direction of 
ferro- and antiferromagnetic order are not independent 
as the ferromagnet already breaks the SU(2) symmetry 
to U(l). In the resulting "spin-flop" phase ferro- and an- 
tiferromagnetic order arc oriented perpendicular to each 
other as has been pointed out in the cold atom context in 
Ref. [ll]. In the SU(3) case, the larger symmetry group 
implies that the staggered order can still point in an ar- 
bitrary direction within the two-component subspace. 



E. Color density wave (CDW) 

While the previous phases can be obtained by both of 
the numerical methods we used, the more general min- 
imization routine found a state with parameters mo = 
3p - 3/2, Ai = A2 > and A3 < 0. This implies that 
the third component is pinned to half filling and occupies 
dominantly a different sublattice than the first two com- 
ponents, leading to a staggered modulation of the total 
density (see Fig.[T|). Due to the doubling of the unit cell, 
the third component is in a band-insulating state. As a 
homogeneous phase, we find that a metallic CDW phase 
is obtained as a minimum in a tiny doping regime close 
to p = 1/2 (see Fig. [3]). Yet we will show below that this 
regime is unstable with respect to phase separation and 
only an insulating CDW state with p = 1/2 is realized. 

At half filling, p = 1/2, the color density wave state 
has been found previously both by fRG on a square lat- 
tice [1^ and by DMFT The energy and staggared 
moments we find at p = l/2,U = 2.5W in the CDW 
phase are consistent with the DMFT calculation [l^ 
within a few percent. 

Note that away from half filling, the color density 
wave, that is, a staggered order with TO3 7^ toi = m2, 
would also imply a uniform polarization of the system, 
P3 7^ Pi = P2- This should be compared to an SU(2) 
antiferromagnet, where staggered order does not induce 
ferromagnetic order. Technically, this difference between 
the SU(3) and SU(2) cases arises because there is no 
symmetry transformation which maps the relevant Gell- 
Mann matrix Ag to — Ag (the two matrices have a different 
spectrum). 

While the fRG for the D = 2 square lattice [l^ seems 
to support that at half filling, the CDW is the ground 
state, in the DMFT study [l^ it is claimed that in the 
SU(3) case these states are degenerate. Within the Gutz- 
willcr calculation we found that at p = 1/2, the density 
wave state has a slightly lower energy than the AFMM 
state, favoring the results from the fRG. However, for 
strong U (even away from half filling) the whole varia- 
tional energy surface Ey = Udtoti^ + 0{W/U)) becomes 
flatter and flatter. Since different states can have the 
same interaction energy oc c?tot — (1/2) + dsa), 

the actual order parameters are determined by the sub- 
leading kinetic energy term. Thus the energy difference 
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between the CDW and AFMM states decreases, and it 
is very difficult to track tlie pliase transition line numer- 
ically for U/W oo. On general grounds, however, it 
is clear that the AFMM can never be the ground state 
at half filling as the Fermi surface of the gapless third 
component shows perfect nesting for a particle-hole sym- 
metric model implying that the state has to be unstable. 

F. Phase diagram 

In Fig. H] wc show which of the states has the low- 
est energy. As all relevant phase transitions are of first 
order, it is, however, clear that the phase diagram has 
to be modified to take phase separation into account. 
Close to a first-order transition one can lower the energy 
by allowing for heterogeneous phases, that is, by mixing 
phases with different densities. While in electronic sys- 
tems macroscopic phase separation is not possible due to 
the presence of long-range Coulomb interactions, it will 
occur for cold atom realizations of the SU(A^) Hubbard 
model. The true phase diagram is obtained by using the 
well-known Maxwell construction for first-order transi- 
tions. For N = 2 the relevant first-order transitions and 
coexistence regions have been discussed in Ref. p^ . 

In Fig. [5] we show the phase diagram after phase sep- 
aration is taken into account. Comparing this to Fig. |4] 
one realizes that both the metallic two component anti- 
fcrromagnct (AF2) and the metallic color density wave 
(CDW) are completely wiped out by phase separation. 
In these regions, the commensurate AF2 insulator (and 
the insulating CDW) coexist with other metallic phases. 
A complex interplay of various phase separated regions 
occurs close to the ferromagnetic transition (see inset of 
Fig. [5]). Note that when two different coexistence regions 
meet, the volume fractions of the phases jump suddenly 
(assuming global thermal equilibrium) when, for exam- 
ple, U is changed. In contrast, when entering a phase 
separated region starting from a uniform phase, the vol- 
ume fractions change always smoothly. 



VI. CONCLUSIONS 

In this work we studied the ground state of the SU(3) 
Hubbard model using a Gutzwillcr calculation. We found 
five phases compatible with two-sublattice symmetry 
breaking on the Bethe lattice. The resulting rather com- 
plex phase diagram is characterized by first-order transi- 
tions and various coexistence regimes. We also compared 
our results to two independent DMFT studies and found 
reasonable quantitative agreement. While we have stud- 
ied the Bethe lattice in the limit of large coordination 
number, z — >■ 00, we expect that the topology of the 
phase diagram will be very similar on a cubic lattice in 
three dimensions. 

It is an interesting question to discuss how the phase 
diagram would manifest itself in a cold atom experiment 




P 



FIG. 5: (color online) Phase diagram based on the Gutzwiller 
calculation. All phases are separated from each other by 
first-order transitions and corresponding coexistence regions 
(shaded). The inset zooms into a region of the phase dia- 
gram where several coexistence regions meet leading to sud- 
den jumps in the volume fractions of the corresponding phases 
as a function of U. The dashed lines represent the stability 
regions of the homogeneous phases shown in Fig. [J] but have 
no direct physical significance. 



assuming that low enough temperatures can be achieved 
to realize all phases. At least two aspects have to be con- 
sidered: First, a parabolic trapping potential holds the 
atomic cloud together in typical experiments. Second, 
due to the underlying SU(3) symmetry, the number of 
atoms of each component is conserved and one can there- 
fore not obtain a state with a finite net polarization start- 
ing from an unpolarized state. This is especially relevant 
as not only the ferromagnetic phase but also the antifer- 
romagnetic phases AF2 and AFMM are characterized by 
a finite net polarization. This implies [H, [s^l that do- 
mains have to form such that the total net polarization 
of the system vanishes. The direct detection of such do- 
mains is possible by taking color-selective phase-contrast 
images of the cloud as long as the domain sizes are larger 
than the spatial resolution jsj. Furthermore, staggered 
order can, for example, be detected by measuring noise 
correlations in a time-of-fiight experiment (35| or more 
directly by Bragg scattering [sl] . While such approaches 
can distinguish, for example, two-sublattice from three- 
sublatticc order directly, multiple domains make it very 
difficult to separate a color density wave from a two- 
sublattice antiferromagnet by measurements of correla- 
tion functions. One can, however, use the fact that the 
color density wave is the only phase with two-sublattice 
structure but no net polarization to distinguish it from 
the other phases. 

To understand the properties of the system in a trap 
it is useful to redraw the phase diagram as a function 
of the chemical potential, fi = dE/3dp, instead of p (as 
shown in Fig. [S]). For a large number of atoms and a 
smooth confining potential V{r) one can locally approx- 
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FIG. 6: (color online) Phase diagram as a function of the 
chemical potential /i — dE/3dp. (Note the logarithmic scal- 
ing.) While the Mott insulator region is rather robust as the 
Mott gap is large ~ U, the gap of the CDW insulator is small 
0{W). All transitions are first order with the exception of 
the transition from AFI to AFMM. The AFI phase can host 
different magnetic structures, see the text. 

imate the inhomogeneous system by a homogeneous one 
with a local chemical potential fi — V{r) (the so-called 
'local density approximation'). From the phase diagram 
shown in Fig. [5] one can therefore directly read off the se- 
quence of phases expected in the presence of a trapping 
potential. 

For large U, the most dominant phase is obviously the 
antiferromagnetic insulator. As it is characterized by a 
large charge gap of order U it is stable for a large range of 
(local) chemical potentials. The other insulating phase, 
the color density wave obtained for p = 1/2 has, in con- 
trast, only a small charge gap which never gets larger 
than ~0.2 W. The difference of the two insulators can be 
understood by considering the limit of vanishing hopping. 
While for p = 1/3 the motion of all particles is locked, 
this is not the case for p = 1/2. Therefore one expects 
that for p ~ 1/2 any gap will be of order W rather than 
U. 

The phase diagram presented here is certainly not com- 
plete. Especially in the AFI phase at p = 1/3, we expect 
that phases we did not discuss appear. For large U/W, 
the insulating phase at low energies is described by the 
SU(3) Heisenberg model. It has been shown [13] that 
the classical ground state is highly degenerate and also 
within the Gutzwillcr approximation a large number of 
ordering patterns have the same energy to leading order 
in W/U. Quantum fluctuations (not described by the 
Gutzwiller wave function) can stabilize certain ordered 
states, and it has indeed been shown by Toth et al. [l^ 
that they favor a three-sublattice ordering relative to the 
two-sublattice antiferromagnetic order discussed by us. 
Higher order terms in W/U may lead to a more com- 
plex magnetic phase diagram within the Mott insulating 
AFI phase, but we expect that this will not change qual- 



itatively the interplay of the AFI phase with the other 
phases at finite doping. For example, a simple argument 
suggests that particle doping favors two-sublattice order- 
ing (the AFMM phase) relative to three-sublattice order- 
ing for large U: When a single 'blue' particle is added to 
a 'red-green' two-sublattice antiferromagnet, it can gain 
the full kinetic energy —W. This is, in contrast, not 
possible when a particle is added to a three-sublattice 
antiferromagnet where all three colors take part in the 
magnetic order. It is also unlikely that subtle quantum 
fluctuations can remove the strong first-order transition 
from the AFI to the PM phase obtained from our varia- 
tional study. 

Overall, the physics of the SU(3) Hubbard model turns 
out to be very different from the well-established SU(2) 
version in many aspects, (i) While in the spin-1/2 Hub- 
bard model the Mott insulating state occurs at half fill- 
ing, where the paramagnetic Fermi surface is perfectly 
nested, the Mott insulator of the SU(3) model at p = 1/3 
shows no nesting in the paramagnetic state, (ii) Due to 
its larger symmetry, the Mott insulating phase can sup- 
port a larger set of competing phases, for example, those 
with three-sublattice structure discussed above, (iii) The 
SU(3) Mott insulator reacts very differently to particle 
doping. Adding more particles does not destroy magnetic 
order as efficiently as in the SU(2) case: two-sublattice 
antiferromagnetism is not frustrated by doping as the 
dynamics of a third species does, to leading order, not 
perturb a magnetic state formed by the first two colors. 
This leads to the stable AFMM phase where two compo- 
nents remain exactly at half filling, forming an orbitally 
selective Mott insulator, (iv) Hole doping, in contrast, 
shows different physics characterized by a strong first- 
order transition to the paramagnetic state, (v) Also, 
all other transitions from the paramagnetic state to or- 
dered states are strongly first order and associated with 
large jumps of the occupations (with the exception of the 
weakly interacting regime very close to half filling), (vi) 
Finally, the physics at half filling, p = 1/2, is fundamen- 
tally different. While a phase with a charge gap is also 
obtained in the SU(3) case, the color density wave does 
not arise from a Mott insulator: its gap remains finite in 
the large U limit. 

For the future, it is an interesting question how more 
exotic phases, like the "chiral spin liquid" [13, HI], ex- 
pected for larger TV, react to doping and whether also 
new types of superconductivity can be realized in repul- 
sive SU(iV) Hubbard models. 
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Appendix A: Cavity theory 

We shall discuss how to calculate local quantities 
within the Gutzwiller approximation using functional in- 
tegrals and the cavity method of dynamical mean field 
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theory. As a concrete example, we shall calculate the 
Gutzwiller expectation value of triple occupancies in the 
limit z oo. Using the definition of the triple occupancy 
operator, we get 



{G\ 



niini2n. 



.3|G) 



(G|G) 



(G|G) 



(Al) 



where yi = E/'c/(^l)'^' Here a creation opera- 

tor is always to the left of the corresponding annihila- 
tion operator Cia , and normal ordering can be performed 
easily. By applying Wick's theorem to the normal or- 
dered expression, we can express both the numerator and 
the denominator as a functional integral. 



(A2) 



We can simply reexponentiate all terms j ^ i in the 
product, however, the term i is missing to complete the 
action S'aux [defined by Eq. ([9])]. To this end we define 
coefficients wj such that 



(A3) 



holds. It can be shown that wj can always be expressed in 
terms of Xj, for example, — wi = (Ag — Af)/A0, 

however, we shall never use these values explicitly. 

The Gutzwiller expectation value of the triple occu- 
pancy is therefore finally expressed as an expectation 
value in the auxiliary theory: 



X^{i)niini2hiz wi{i)nu 



(A4) 



Now we integrate over all Grassmann variables except for 
the site i, which we call the cavity. Formally, this gives 



ni 



(A5) 

however, the calculation of the cavity action 5cav[*] is not 
possible in general. But in the limit of z oo. one finds 
that the contributions from other sites just renormalize 
the quadratic terms [111, and therefore the cavity action 
is 



I 



(A6) 



where the cavity bare propagator D'^ has to be deter- 
mined self-consistently from the condition that any local 
quantity defined at site i in the auxiliary theory has to 
have the same expectation value as in the cavity theory. 
We shall postpone the calculation of T)'^ for the moment. 

The main advantage of the cavity method is that we 
can calculate any finite dimensional Grassmann integrals 
in the cavity explicitly, that is, solve the "impurity prob- 
lem" exactly. In contrast to dynamical mean-field theory 
this is always possible since the action is static. The 
cavity partition function is simply given by 



1 



(A7) 



where we introduced the notation 



A, = A2(1-I?0)(1-I?0)(1-P3"), 

A, = A?2?°(l-2?2")(l-2?°),... 

A-, = Af(l-2?0)2?2"2?3",... 

At ^ XlV°V°V°. (A8) 

Finally, we can perform the Grassmann integrals in 
Eq. (|A5P and find that the triple occupancy is simply 



(A9) 



To derive this expression, we do not need the explicit 
expressions for uj and wj in terms of the Gutzwiller pa- 
rameters A/, just the fact that both the exponential rep- 
resentation e^^i^'^' and the "inverse" '^jWihi exist. 

Repeating analogous steps for the double and single 
occupancies lead to expressions, for example. 



{G\n,2n^3\G) Ai{i) + At{i) 



(G|G) 



E/^/« 



(AlO) 



and 



{G\hn\G) _ A,{i) + A2{t) + A-^ii) + At{i) 



(G|G) EiAiit) 
Let us now fix A0 and A^ such that 



(All) 



A1 + A2 + A3 + At 
A2 + A1 + A3 + At 
A3 + Ai + A2 + At = P^, 



(A12) 



and discuss the implications. 

First, we shall determine 2?° from the self-consistency 
relation ("DMFT equation") for the local propagator. 
This means that 



Gi; 



V\i] 



(A13) 
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The left-hand side can be calculated using Dyson's equa- 
tion on the lattice, 



(A14) 



where E is the proper self-energy matrix in the auxil- 
iary field theory, which becomes site diagonal in the limit 
2; — )■ oo [m . The right-hand side of Eq. (|A13|) can be cal- 
culated analytically and one finds that 



If we use the conditions in Eq. (jA12p . we simply get 



(A15) 



V — V 



(A16) 



But this means that the proper self energy of the cavity 
theory vanishes. Furthermore, in the limit z — > 00 the 
proper self energy in the cavity is equivalent to the self 
energy on the lattice [2^, and therefore 



(A17) 



From Eqs. ( \^ ~ (|llT|) we get 

= = <, (A18) 

= A2(z)(l-On02n03-f A?(zXn°2n°3,(Al9) 
U = \\(S)n\^n\^T{i^. (A20) 



We conclude that fixing the single-occupancy param- 
eters Xa in such a way that the physical and the bare 
densities coincide, implies that the proper self energy in 
the auxiliary field theory vanishes for z — >■ 00. 



Also, Eqns. (dU]), (HH, (|A19|) . and (|A20|) can be used 
to replace the Gutzwiller variational parameters by the 
physical occupancies. The relations take the form of 
"mass laws" [2J|, for example, 



\\(^) 



K0 

1 - 



HI 



HZ 



rP 
Pil 



H2 



Pil da - U 



Pit _ u 

PU <"°2"°3 ' 



(A21) 
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